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Spacecraft in orbit near libration point Li in the Sun-Earth system are 
excellent platforms for research concerning solar effects on the terrestrial 
environment. One spacecraft mission launched in 1978 used an Li orbit for 
nearly 4 years, and future Li orbital missions are also being planned. Orbit 
determination and station-keeping are, however, required for these orbits. In 
particular, orbit determination error analysis may be used to compute the 
state uncertainty after a predetermined tracking period; the predicted state 
uncertainty levels then will impact the control costs computed in 
station-keeping simulations. Error sources, such as solar radiation pressure 
and planetary mass uncertainties, are also incorporated. For future missions, 
there may be some flexibility in the type and size of the spacecraft’s nominal 
trajectory, but different orbits may produce varying error analysis and 
station— keeping results. The nominal path, for instance, can be (nearly) 
periodic or distinctly quasi-periodic . A periodic "halo" orbit may be 
constructed to be significantly larger than a quasi-periodic "Lissajous" path; 
both may meet mission requirements, but perhaps the required control costs for 
these orbits are provably different. Also for this spacecraft tracking and 
control simulation problem, experimental design methods can be used to 
determine the most significant uncertainties. That is, these methods can 
determine the error sources in the tracking and control problem that most 
impact the control cost (output); it also produces an equation that gives the 
approximate functional relationship between the error inputs and the output. 
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INTRODUCTION 

In one formulation of the problem of three bodies, when the mass of one of 
the bodies is sufficiently small (infinitesimal) so that it does not affect the 
motion of the other two bodies (primaries), the "restricted three-body problem" 
results. Five libration (Lagrange) points can be found as particular solutions of 
the equations of motion governing the path of the infinitesimal mass moving within 
the gravitational fields of the primaries. These equilibrium points are defined 
relative to a coordinate system rotating with the primaries. One Lagrange point, 
Li, is located between the primaries and is the libration point of interest here. 

Three-dimensional, periodic and quasi-periodic orbits are currently being 
studied for upcoming missions. Periodic "halo" orbits in the vicinity of 
libration points have been studied since the late 1960s. Robert Farquhar coined 
the term "halo" to describe a three-dimensional, periodic orbit near a libration 
point on the far side of the Moon in the Earth-Moon system. These orbits were 
designed to be large enough so that the spacecraft would be constantly in view of 
the Earth and thus would appear as a halo around the Moon. Alternatively, the 
variations in size and shape that a quasi-periodic orbit can exhibit may add 
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valuable flexibility for mission planning. This type of bounded, 
three-dimensional libration point trajectory is called a "Lissajous" orbit since 
specific planar projections of these quasi-periodic trajectories may look like a 
special type of Lissajous curve. 

2 

Howell and Pernicka have developed a numerical technique for determination 
of three-dimensional, bounded libration point trajectories of arbitrary size and 
duration. Their numerical algorithm uses an analytic solution as a first 
approximation and then constructs a trajectory continuous in position and 
velocity. Their method is used in this study to define nominal paths in the 
Sun-Earth+Moon problem. The notation "Earth+Moon" means that the Earth plus the 
Moon are treated as one body with mass center at the Earth-Moon barycenter. The 
numerical data is then curve fit using a cubic spline routine, although the use of 
other curve fit methods is possible. The assumed dynamic model is the elliptic 
restricted three-body problem (ER 3 BP), where the primaries move on known elliptic 
paths. The force model used here includes solar radiation pressure 4 , the 
gravitational attractions of the Sun and the Earth+Moon barycenter, and the 
centrifugal force associated with rotation of the system. 

The forces affecting the spacecraft orbit have differing levels of 
uncertainty, and, unfortunately, the spacecraft will drift from the nominal path. 
Both range and range-rate tracking also include inaccuracy in measurement. The 
accumulated error in the spacecraft’s position and velocity relative to the 
nominal path after a predetermined period of tracking can be computed. This 
error, or uncertainty, in the spacecraft state is measured through simulations, 
commonly referred to as orbit determination error analysis, and is presented as 
a vector of standard deviations of the states. In this work, the state vector 
includes three position and three velocity states. The state uncertainty computed 
in the error analysis can then be input to a station-keeping algorithm that 
computes control manuevers to return the spacecraft to the vicinity of the nominal 
path. The algorithms incorporate certain minimal constraints for time between 
manuevers, control magnitude, and distance from the nominal path before a control 
manuever is input. For these algorithms, variations in orbital shapes and sizes 
may have some effect on the station-keeping costs. 


Coordinate Systems 

The coordinate systems used 
in this analysis have a common 
origin at the primaries’ center 
of mass. Primaries with masses 
mi and m2 such that mi £ m 2 are 
assumed here. The infinitesimal 
mass is denoted as m 3 .. These 
masses (mi, m2, m3) correspond to 
particles situated at points Pi, 
P2, and P3 , respectively. The 
barycenter is denoted as "B," 
and the resulting arrangement is 
shown in Fig. 1 . The rotating 
coordinate system is defined as 
XrYrZr, and the inertial system 
is identified as X1Y1Z1. 


BACKGROUND 



Fig. 1 Coordinate Systems Used 
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Note that both coordinate systems are right-handed, and the X and Y axes for both 
systems are in the plane of motion of the primaries. The rotating Xr axis is 
defined along the line that joins the primaries and is directed from the larger 
toward the smaller primary. 


Equations of Motion 

The equations of motion for m3 (the spacecraft) relative to B as observed in 
the inertial reference frame are now formulated. The sum of the forces on m3 
resulting from both the gravity fields of masses mi (the Sun) and m 2 (the 
Earth-Moon barycenter) and from the solar radiation pressure can be used to 
produce the following second-order vector differential equation: 

p" = - G d - G (-^-) ? + (-^f-) d. (1) 

d 3 r 3 d 3 

The overbar denotes a vector, and primes indicate differentiation with respect to 
dimensional time. All quantities are dimensional, as appropriate, and ^ the 
quantity "G" is the universal gravitiational constant. The scalars d and r in 
Eq. (1) denote the magnitudes of vectors d and r, respectively, depicted in Fig. 
1. The dimensionless scalar "k" is the solar reflectivity constant, and S is 
the solar radiation pressure constant . The position vector p is written in 
rotating components as 

p = x Xr + y Yr + z Zr (2) 


where Xr,Yr,Zr are unit vectors. The kinematic expression for p" is: 

-« = (x"-e"y-2e'y'-e ,2 x) Xr + (y"+ 0 "x+ 20 , x'- 0 ,2 y) Yr + z" zr. 


(3) 


Three scaled equations of motion for P3 can be derived using the following 
definitions: the sum of the primary masses is one mass unit, the mean distance 

between the primaries is one distance unit, and the universal gravitational 
constant is equal to one unit by proper choice of characteristic time. The 
equations of motion can then be simplified and scaled by also introducing the 
nondimensional mass ratio p, "psuedo-potential" U, and the scaled solar radiation 
constant s: 


H = 


m 2 

mi + m2 


u = * 

a 


1 *2,2 2 , k S 

+ -—e(x + y ) - — 


(4) 

(5) 


where the dot denotes the derivative with respect to characteristic time. Then 
the vector magnitudes, "d" and "r," are written in terms of scaled quantities as: 


r . _ . 2 2 2 , 1/2 
d = [ (x + p R) + y + Z ] 

,2 2 2,1/2 
r=[(x-R+pR) + y +z] 


(6) 

(7) 


The three second-order differential equations that result can be written in 
terms of characteristic (scaled) quantities as 

x - 2 0 y = 4^- + 0 y = Ux + 0 y, (8) 


y + 2 0 x = 


9x 

au 

ay 

au 

az 


- 0 X = Uy - 0 X, 
= Uz. 


(9) 

( 10 ) 


These three equations can then be used to propagate the spacecraft state forward 
in both the error analysis and station-keeping simulations. 


397 



Reference Paths Generated for This Work 


In the ER3BP , precisely periodic halo orbits exist, but nearly periodic 
orbits are more practical and likely to be used in mission planning. Therefore, 
the goal here will be to compare results for quasi-periodic Lissajous and nearly 
periodic "halo-type" orbits. Fig. 2 depicts one orthographic view of the 


Lissajous and halo-type orbits used here 
(approximately four times) larger in both 



y in kilometers xl05 
Fig. 2 Orthographic Depictions 

Curve Fitting the Nominal Path 


The halo-type orbit is significantly 
the X and Y excursions from Li . 



y in kilometers xl05 
of the Reference Trajectories 


A numerical integration method developed by Howell and Pernicka 2 is used to 
generate a set of reference points at specified times for both position (three 
states) and velocity (three states), relative to the libration point of interest. 
The method computes numerical data in a reference frame that is centered at the 
libration point (in this case Li ) and rotates with the primaries. However, the 
state estimation techniques and station-keeping algorithms used in this work 
require access to a continuous nominal path of acceptable accuracy. In one study, 
Pernicka found that station-keeping costs for a libration point orbit were 
sensitive to the accuracy of the curve fit. A cubic spline interpolation routine 
was selected to model the reference trajectory here; the results of using other 
methods are summarized in the station-keeping section of this effort. 

Examples of Experimental Design (DOE) Methods 

DOE methods are used to purposefully change the most important inputs to a 
process in order to analyze the output. The inputs are coded and alternately set 
at predetermined values for each experimental run so that the design is 
orthogonal; the relative contribution of each input can thus be judged 
independently. The output of interest may be the mean response and its variation, 
with the ultimate goal being to hit an output target value and minimize the output 
variability. However, the results from the set of experimental runs also 
determine the estimated function that relates the inputs to the output (s). 
Experimental design methods are also used to reduce the required number of runs or 
screen out relatively unimportant input variables. When only three inputs at two 
different input levels are considered, a two-level, "full-factorial" design 
consisting of every possible combination of input factors would require 2 3 = 8 
total runs. This design allows the experimenter to obtain the full model with all 
possible interactions. However, if 7 inputs in a 2— level design were used, 

2 — 128 individual runs would be necessary. These 128 runs may be expensive in 
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terms of both time and money. As a result, fewer runs may certainly be desired. 
A full-factorial two-level design will include all input variables and their 
interactions in the output equation model, but this may not be required. If some 
interactions are known to be relatively unimportant, a fractional-factorial 
design consisting of a fraction of the number of runs required in a full-factorial 
can be constructed. For instance, for 3 inputs, a full-factorial would 

necessitate 2 3 = 8 runs, while one type of fractional-factorial would require as 
few as 4 runs to determine the significance of main effects (modeling no 
interactions) . 

In a recent work by Garrett 8 , several simple, yet educational, examples of 
the use of DOE are described. A similar example is included here: it is assumed 

that the area of a rectangle can be measured precisely, but the functional 

relationship between area and the length and width is not known. This example is 

truly hypothetical, but it can be used to illustrate simple DOE computations. The 
"design space" (where the computed model can be considered a good approximation to 
the true system) is defined by 1 ^ width ^ 2 meters and 1 ^ length ^ 3 meters. 
Here, M w" is used for width and "l" for length. Runs are accomplished at the 
extreme values of the input variables, with w = 1 or 2 meters and l—l or 3 

meters; however, first these measurements are generally coded. The data is coded 
by using the averages of both measurements and their ranges (highest value minus 
lowest). With R(l) = range of l t R(w) = range of w, w = average of the w 
extremes, and I = average of the l extremes, the coded settings are Wc and lc : 


wc=2 




( 11 ) 


When coded, the extreme values become +1 and -1 for each input, and these values 
are more simply denoted as "+" and respectively. A balanced design with 4 
runs then yields a design matrix of 


RUN Wc lc (Wc)(£c) 

1 + 

2 - + 

3 + - 

4 + + + 

The experiment is conducted using these high and low settings, and the measured 
areas of the rectangle (outputs) are 1, 3, 2, and ? 6 square meters for runs 1 

through 4, respectively. Schmidt and Launsby discuss interesting hand 
computational methods to determine the output equation; however, simple least 
squares methods also provide identical results. The prediction equation for the 
output is assumed to be 

a = bo + bi Wc + b2 lc + b3(wc) (£c) (12) 


where a = estimated area and the coefficients are computed using a least squares 
method with 


This method yields 


■ 1 

-1 

-1 

1 


■ i 

1 

-1 

1 

-1 

and a = 

3 

1 

1 

-1 

-1 

2 

1 

1 

1 

1 


6 


b = (C T C) _1 C T a = [bo bi b2 b3] T . 

a = 3 + Wc + 1 . 5 tc + .5 (wc) ((c) . 


in 


(13) 

(14) 


A similar method could be used to derive a prediction equation for the ^variance 
(or the natural logarithm of the standard deviation) of the output . The 
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resulting model in Eq. (14) is normally checked by completing test (confirmation) 
runs at extreme values and at the midpoint of the design space. For this example, 
that could mean using values of +1 for both wc and tc for an extreme run or 0 for 
both inputs for the midpoint run. Suppose a confirmation run was conducted at the 
extreme values and the resulting output was a = 6 meters. Using Eq. (14), the 
predicted output is 

a = 3 + 1 + 1.5(1) + .5(13(1) =6 meters. 

For a confirmation run at the midpoints, the measured answer is a = 3 meters. 
Using Eq. (14), the predicted output equation is 

a = 3 + 0 + 1.5(0) + . 5(0) ( 0 ) = 3 meters. 

Hence, the confirmation runs verify the model; a significant disagreement would 
require further investigation. (In fact, this Is the exact functional model — it’s 
just coded. ) When noise in the system exists, statistical tests are used to 
test confirmation. The coded Eq. (14) can now be converted to use uncoded inputs 
by using Eq. (11 ) : 

= (w )U). 

This example was simplified because we obviously knew the actual output 
equation. In manufacturing or engineering problems, the relationship between 
inputs and outputs is only generally known, and DOE can be used to gain problem 
insight. In the next section, the orbit determination error analysis methods used 
in this effort are summarized. The following section describes the 
station-keeping methods derived for this work and summarizes the control-cost 
comparisons of halo and Lissajous orbits. Finally, modeling the inputs of the 
station-keeping routine in an experimental design is presented. 

ORBIT DETERMINATION ERROR ANALYSIS 

Complete, exact knowledge of the state of a spacecraft in orbit is generally 
not possible. Available measurements are usually some function of the state 
variables and are not precise. For instance, a spacecraft in a libration point 
trajectory in the Sun-Earth system may be tracked using range and range-rate 
measurements containing random errors. The spacecraft may be affected by forces 
inadequately represented in the dynamic model, and model parameters may be 
uncertain. By definition, the linearized system of equations used to model the 
nonlinear system is a further approximation. These sources of error make 
knowledge of the spacecraft state uncertain. Computation of the most likely 
current state of the spacecraft in the presence of measurement and model 
uncertainty is the focus of orbit determination. 

Error analysis involves an investigation of the impact of various error 
sources on orbit determination. The outputs of this error analysis are the 
standard deviations of the states. These outputs could then be used to predict 
how an improvement in measurement accuracy, for instance, would lessen state 
uncertainty. One benefit of more accurate knowledge of the state might be a 
reduction in station-keeping costs. A mathematical procedure can be developed to 
combine all information about the spacecraft state, filtering this observed data 
based on the varying degrees of uncertainty, to obtain a "best estimate" of the 
state and an estimate of the resulting state variable uncertainties. 

The measurement and dynamic models are first summarized, three error analysis 
methods are briefly discussed, and then results are summarized. The three error 
analysis methods used here are the^ Kalman filter, batch weighted least squares , 
and consider covariance analysis ' . Each technique computes a covariance matrix 


A 

a 


= 3 + (2) (w - 1.5) + (1.5) (2) 




+ ( . 5) (2) (w-1 . 5) (2) 
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at a specified epoch, and the positive square roots of the diagonal entries are 
indicators of state uncertainty levels. 

Measurement and Dynamic Models 

The measurement and dynamic models used in the filter derivations are 

Measurements: Zk = Mk Xk + Vk, (15) 

Dynamic Model: Xk+i = $(tk+i,tk)Xk = $(k+l,k)Xk, (16) 

where Zk is the measurement residual vector at time step k; Xk is the residual 
state vector at time step k; Mk is the measurement matrix that is linearized about 
the nominal path; $(k+l,k) is the state transition matrix at time step k+1 
relative to time step k; and Vk is the measurement noise vector with_assumed 
statistics E(Vk)=0 and E(VkVk)=Vk, where "E” is the expectation operator, 0 is the 
zero vector, and Vk is the measurement noise covariance matrix. Range and 
range-rate measurements are assumed; the matrix, M, is then a time-varying matrix 
of dimension 2x6, evaluated along the nominal path. 

Error Analysis Methods Used 

Early work in this area was designed to compare the error levels obtained 
here to those found in other works and to determine error levels for use in 
follow-on station-keeping simulations. Three methods of orbit determination error 
analysis (using covariance analysis) were formulated: Kalman filter, batch 

weighted least squares, and consider covariance analysis. The results of Kalman 
and batch weighted least squares filters were virtually identical, as expected, 
but nonetheless helped to confirm the analysis. Both methods were formulated to 
compute state uncertainty after a predetermined number of tracking updates, 
simulating range and range-rate measurements with associated error statistics. 
Consider covariance analysis also uses a batch weighted least squares formulation 
but includes parameter uncertainty. Model parameters that were initially 
considered uncertain in this work were the planetary masses (through the 
dimensionless mass parameter p), the locations of the tracking stations, and the 
solar reflectivity constant. In general, at the epoch of interest, the state 
uncertainty is considered the consequence of the accumulated uncertainties in the 
model, the parameters of interest, and the measurements 

Orbit Determination Error Analysis Results 

A survey of input error levels used in similar error analysis studies serves 
as a valuable introduction. The values of these uncertainties may be used to 
compute diagonal entries of input covariance matrices for an error analysis, or, 
alternatively, may be used as error sources in a station-keeping simulation. 
Table 1 lists the input error levels assumed in several error analysis studies. 
The errors are denoted by the symbols generally used in the derivation sections of 
this work. The solar reflectivity constant is k; the tracking site location 
uncertainty is S and is input as an equal uncertainty level for each of the site 
coordinates xs, ys, zs; range tracking is R; range-rate tracking is RR; and the 
uncertain mass parameters are pe for Earth, ps for the Sun, and pm for the Moon. 
The last column contains the uncertainty in dimensionless mass parameter p that 
would be "equivalent" to the errors listed for the individual mass parameters. 
(Recall that p = (pe + pm)/(ps + pe + pm) for the three-body system of interest in 
this work.) The approximate value of <r( p) (standard deviation of p) is calculated 
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from extensive (10,000) Monte Carlo trials for each of these studies. An entry in 
Table 1 of " — " means the particular study did not indicate if an uncertainty of 
this type was used. 

Table 1 

SURVEY OF ERROR ANALYSIS INPUT ERRORS 




ONE STANDARD DEVIATION ERRORS 



STUDY 

k 

S 

(km) 

R 

(km) 

RR 

(m/sec) 

M® Ms 

( <r km 3 /sec 2 

fim 

) 

<r(p) 
(x p) 

Mistretta 13 

15% 

— 

.010 



3.08 

xlO 6 

.0726 

2.335 

xlO -5 

Joyce 1 4 

10% 

.002 

.015 


.3986 

1.327 

xlO 4 

.0490 

1.411 

xlO" 7 

r> r 15,17 

Efron 

Rodriguez- 

10% 

.002 

.015 


.3986 


.0490 

1.411 

xlO’ 7 

Canabal 16 

— 

.010 

.015 


— 

— 

— 

— 

Longuski 17 

13% 

. 0003 

.010 


.4903 

4030.7 


1.231 

xlO' 7 

This Work 

13% 

.010 

.015 


.3986 

1.327 

4 

xlO 

.0490 

1.411 

xl0~ ? 


The error analysis conducted here assumes a 20-day tracking arc with 3 passes 
per day from 3 tracking sites. These assumptions closely match those of Joyce 14 . 
Using this tracking schedule and the R and RR measurement errors listed for this 
work in Table 1, the Kalman filter produces error levels presented in Table 2. 

Table 2 

RESULTS USING A KALMAN FILTER FOR ERROR ANALYSIS 
I ONE STANDARD DEVIATION LEVELS OF STATE ERRORS ] 

x (km) y (km) z (km) x (mm/sec) y (mm/sec) z (mm/sec) 


.550 1.600 4.450 .430 .775 2.250 


The error levels listed in Table 2 are a result of a covariance analysis for 
the halo-type nominal path. The magnitudes of the error levels listed in Table 2 
are, in fact, quite small; when additional error sources, such as mass parameter 
and station location uncertainties, are included in a consider covariance 
analysis, the resulting state error levels increase. The results in Table 3 are 
from a consider covariance analysis incorporating R and RR tracking, station 
location, and mass parameter uncertainties at the levels listed in Table 1 for 
this work. 
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Table 3 

ERROR LEVELS PRODUCED FROM CONSIDER COVARIANCE ANALYSIS 



One Standard 

Deviation Levels 

Coordinate 

Halo-Type Orbit 

Lissajous Orbit 

x (km) 

1.46 

1.25 

y (km) 

2.64 

3.35 

z (km) 

4.81 

3. 19 

x (mm/sec) 

1.40 

1.25 

y (mm/sec) 

1.85 

1.41 

z (mm/sec) 

2.49 

2.51 


It certainly may be of great interest to compare the error levels found in this 
effort with the results of other investigations involving spacecraft in halo (or 
halo-type) orbits near the interior Sun-Earth libration point. Table 4 lists the 
results of four studies that do not include solar reflectivity as an error source 
and have small differences in the nominal paths and force models. 


Table 4 

COMPARISON OF ERROR ANALYSIS RESULTS FROM SEVERAL SOURCES 



One Standard 

Deviation 

Error Levels 


Coordinate 

Rodriquez-Canabal 1 6 

,18 

Simo 

,19 

Simo 

This Work 

x (km) 

2.7 

1.5 

1.73 

1.46 

y (km) 

3.9 

2.5 

2.24 

2.64 

z (km) 

3.4 

15.0 

5.48 

4.81 

x (mm/sec) 

2.4 

1.0 

1.41 

1.40 

y (mm/sec) 

3.5 

1.0 

1.41 

1.85 

z (mm/sec) 

1.3 

3.0 

2.45 

2.49 


The differences in error levels listed in Tables 3 and 4 may not be 
statistically significant; that is, station-keeping costs, determined through 
simulations using these error levels, may or may not differ statistically 

The results using one derived control scheme and the data in Table 3 are 

summarized in the next section. 

STATION-KEEPING SIMULATIONS 

For a coll inear libration point orbit, a small deviation from the(unstable) 
nominal trajectory can lead to rather large drift from the path in a short time. 
In effect, a station-keeping algorithm must combat both the current drift from the 
path in addition to the exponential increase in the drift that is expected if no 
correction is implemented. Any delay in the control actuation may allow the drift 
to increase and thus compound the station-keeping problem. The goal of the 

station-keeping routine is then to keep the spacecraft "close enough" to the 

reference trajectory. The allowable deviations may depend on the simulation 
experience with a given control algorithm and on mission constraints, including 
the propellant cost that can be tolerated and the minimal time between control 


403 




inputs. When the spacecraft is "near" the nominal trajectory, it is reasonable to 
model the deviations from the reference path using a linear analysis. 

Derivation of Method 


For the linear control scheme developed here, the state transition matrix is 
partioned into four 3x3 submatrices as 


*(tk,ta) = 


Ak0 

Bk0 

Ck0 

Dk0 


fAk Bk 


Ck Dk 


(17) 


A Av input (a 3x1 vector), with magnitude denoted as Av, is assumed to be added at 
a time to. The Av (delta-velocity) is added to the initial velocity states in the 
numerical integration routine in order to change the deviation of the spacecraft 
from the nominal path at some future time. In this derivation, pk is the position 
deviation (a 3x1 vector) and vk is the velocity deviation (a 3x1 vector) of the 
spacecraft from the nominal path at time tk, with k = 1, 2, 3 and 4. If vo is the 
residual velocity (a 3x1 vector) and pe is the residual position (a 3x1 vector) 
relative_to the nominal path at time to, then a Av input at to could be used to 
predict pk for k = 1, 2, 3 and 4. For instance, when the initial position Xo 
includes an initial velocity perturbation vo, a delta velocity Av, and an initial 
position perturbation po, the state propagation equation results in 



Pk 


P0 

Xk = 

Vk 

= 4>(tk,t0) X0 = $(tk,t0) 

V0+Av 


The cost function used to derive this control scheme is 


J (Av) — Av Q Av + pi R pi + vl Rv vi + p2 S p2 + V2 Sv V2 

+ p3 T P3 + V3 Tv V3 + pi U P4 + V4 Uv V4, (19) 

where Q is a positive definite weighting matrix and R, R v , S, Sv, T, Tv, U, and Uv 
are positive semidefinite weighting matrices. The cost function can be written in 
terms of Av by using substitutions for pk and vk^ with k = 1, 2, 3, and 4, derived 
from Eqs. (17) and (18). The minimum is then Av = 

-[Q+BlRBi + B 2 SB 2 + B 3 TB 3 + B4UB4 + DiRvDi + D2 SvD 2 + D3TvD3 + D 4 UVD 4]" 1 

x[(BlRBl + B 2 SB 2 + B 3 TB 3 + B 4 TB 4 + DiRvDi + D2SvD2 + D3TvD3 + D4UvD4 )v0+ 

(B 1 RA 1 + B 2 SA 2 + B 3 TA 3 + B 4 UA 4 + DiRvCl + D2 SvC2 + D3TvC 3 + D4UvC4)po]. 

A simpler version of this controller can be used by setting, for instance, the 
weighting matrices U and Uv equal to the 3x3 zero matrix. This modified 
controller is the one used in the following section. 


Comparison of Halo-Type and Lissajous Orbits 

The cost of maintaining the spacecraft in orbit for 2 years is selected as 
the comparison value. For each simulation run, tracking updates, with assumed 
srror levels listed in Table 3, are input every 20 days. Solar radiation pressure 
uncertainty is also input as an error source with magnitude listed in Table 1. 
The errors are modeled as zero-mean Gaussian random variables. Each simulation of 
the station-keeping algorithm will be a random trial with the random variable of 
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interest being the total magnitude of the station-keeping costs (AVt) for the 
2-year simulation. A sequence of 30 Monte Carlo station— keeping simulations 
produces a random sample of 30 random variables. Sample statistics, such as means 
and standard deviations, can then be calculated, and statistical tests can be 
conducted to compare the mean control costs for halo-type and Lissajous 
orbits 5 ’ 12 . Table 5 contains the results of one set of simulations using 30 Monte 
Carlo trials for each type of orbit. 

Table 5 

COMPARISON OF STATION-KEEPING COSTS 


Lissajous Orbit 

Halo-Type Orbit 

Avg Avt Std Dev Range 

Avg Avt 

Std Dev 

Range 

(m/s) (m/s) (m/s) 

(m/s) 

(m/s) 

(m/s) 

.8450 . 1603 .57 - 1.15 

.8124 

.1233 

.62 - 1.08 


Statistical hypothesis tests conclude that the 2-year mean control cost, using the 
two nominal paths and this particular controller, are equal. The conclusion of 
equal station-keeping costs for all nominal paths near this libration point and 
any control scheme cannot be drawn from this work. 

Comparison of Station-Keeping Costs for Different Curve Fitting Options 

Various curve fitting methods have been developed to model the nominal paths. 
While cubic splines are used here, least squares curve fits for a trigonometric 
series and linear interpolation routines have also been tested. The data in Table 
6 summarizes efforts to date. The curve fits are indexed by the number of terms 
included in the Fourier series. The cubic spline and linear interpolation schemes 
are indexed by the time between points. 


Table 6 

COMPARISON OF CONTROL COSTS BY CURVE FITTING TECHNIQUES 


Cubic Spline 

Average 2-Year Cost (meters/sec) 

Days between points =3, 6, 9 

1.234, 1.801, 10.324 

[Fourier Series ij 

I Terms Used =28, 91, 121, 161 

9.577, 8.147, 1.419, 1.414 

1 Linear Interpolation 

Days between points = .5, 1, 2, 6 

1.290, 1.307, 1.333, 38.604 


EXPERIMENTAL DESIGN RESULTS 

The process of interest here is station-keeping for a 2-year halo-type 
Lagrange point orbit in the Sun-Earth+Moon elliptic restricted three-body problem. 
The input variables include tracking errors (track), solar radiation pressure 
(SRP) and mass ratio (mass) uncertainties, orbit injection errors (inject), and 
thruster (thrust) errors. The outputs of interest are the 2-year control cost 
(AVt) and its variance. Other inputs could be considered, and additional outputs, 
such as the number of AV inputs required or the average separation time between 
control inputs, could also be evaluated in future efforts. The relationship of 
the inputs, the process, and the outputs is depicted in Fig. 3. 
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INPUTS 


PROCESS OF INTEREST 


OUTPUTS 


TRACKING ERRORS (TRACK) 

SOLAR RADIATION PRESSURE UNCERTAINTY 

7 

(SRP) 

MASS PARAMETER UNCERTAINTY (MASS) 

> 

ORBIT INJECTION ERRORS (INJECT) 

> 

THRUST INPUT ERRORS (THRUST) 

} 

> 


STATION-KEEPING 


2-YEAR 
CONTROL 
COST (AVt) 

AND ITS * 

VARIANCE (S 2 (AVt) ) 


Fig. 3 Process of Interest 


For this analysis, a fractional factorial two-way design was selected in 
order to limit the total number of runs. A fractional factorial 2 5 ” 1 design 7 
allows use of only 16 runs to pick out contributions of the 5 main inputs and 10 
two-way interactions. The design matrix, with only the main effects listed, is 
depicted in Fig. 4. 

RUN INPUT VARIABLES 



a 

b 

c 

d 

e=abcd 


TRACK 

THRUST 

SOLAR 

HASS 

I NJECT 

1 

- 

- 

- 

- 

+ 

2 

- 

- 

- 

+ 

- 

3 

- 

“ 

+ 

- 

- 

4 

- 

- 

+ 

+ 

+ 

5 

- 

+ 

- 

- 

- 

6 

- 

+ 

- 

+ 

+ 

7 

- 

+ 

+ 

- 

+ 

8 

- 

+ 

+ 

+ 

- 

9 

+ 

- 

- 

- 

- 

10 

+ 

- 

- 

+ 

+ 

11 

+ 

- 

+ 

- 

+ 

12 

+ 

- 

+ 

+ 

- 

13 

+ 

+ 

- 

- 

+ 

14 

+ 

+ 

- 

+ 

- 

15 

+ 

+ 

+ 

- 

- 

16 

+ 

+ 

+ 

+ 

+ 



Fig. 4 

Design 

Matrix 



A full factorial would enable analysis of 5 main effects, 10 two-way 
interactions, 10 three-way interactions, 5 four-way interactions, and 1 five-way 
interactions. Generally, interactions above two-way are not significant 
contributors to a model . The modeled interactions not depicted in Fig. 4 are ab, 
ac, ad, ae^ be, bd, be, cd, ce, and de. Note also that the main effect "inject" 
is aliased with the abed four-way interaction. The full factorial two-level 
design would allow analysis of all possible interactions that could affect the 
output. It would not allow curvature analysis (quadratic effects), but these 
could be analyzed using a sequential central composite design approach 7 . 
Investigation of quadratic effects would be necessary only if confirmation runs 
indicate poor agreement at the midpoint of the design space. The design space is 
determined by the extreme values selected for each input. That is, the low and 
high settings for each input determine the region over which the approximate 
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output equations are defined. A large range for an input will also have a great 
bearing over whether it will be found significant. For this design, the low and 
high settings for each input are depicted in Table 7. 


Table 7 

LISTING OF DESIGN INPUTS AND SETTINGS 


Input 


Settings 



Low (-) 

High (+) 

Track 

X 


1 KM 

3 KM 

y 


2 KM 

10 KM 

z 


4 KM 

15 KM 

x dot 


.0010 m/sec 

.008 m/sec 

y dot 


.0015 m/sec 

.010 m/sec 

z dot 


. 0020 m/sec 

.015 m/sec 

Solar 


2 . 5% 

15% 

Mass 


1 . 231xl0 -7 

5. 000x10" 6 

Inject 

Each position 

coordinate 

1.5 KM 

50 KM 

Each velocity 

coordinate 

.001 m/sec 

.05 m/sec 

Thrust (each direction) 

2.5% 

10% 


These input settings are representative of those used in other orbit 
determination error analysis and station-keeping studies. (See Tables 1 and 4. ) 
The resulting output equations for the predicted AV (here denoted as AV) and 
natural logarithm of the output variance ( denoted as ln(S)) are 

A^ = 2.9348 + 1.3627 Track + .2953 Thrust + .0263 Solar + .0028 Mass 

- .0904 Track-Thrust -.0605 Track-Solar -.0215 Thrust-Solar 

- .0103 Track-Mass - .0065 Thrust-Mass + .0395 Solar-Mass + 

- .0090 Solar-Inject + .0551 Thrust-Inject + .0094 Mass-Inject 

- .0299 Track-Inject + .05 Inject, (20) 

ln(§) = 0.2027 + 0.5948 Track + .2334 Thrust + .1221 Solar + .0567 Mass 

- .1614 Track-Thrust -.0959 Track-Solar + .0093 Thrust-Solar 
+ .0517 Track-Mass + .0074 Thrust-Mass - .0297 Solar-Mass 

- .0339 Solar-Inject - .0394 Thrust-Inject 

+ .0238 Mass-Inject - .0127 Track-Inject + .0292 Inject. (21) 

Additional experimental runs showed that the output model confirmed at the design 
midpoint and at both extremes. Often, this sort of model is used to determine 
optimal input settings: in order to minimize both AV and ln(S) in Eqs. (20) and 

(21), all inputs should be set at the minimum settings. However, a more realistic 
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use of these equations is for sensitivity analysis: the size of the coefficient 
of each input is a measure of that variable’s influence on the output. These 
results show that tracking and thrust input errors are responsible for a large 
portion of the control cost. By reducing these two errors to their minimum, 
nominal savings on the order of 50% are predicted. 

CONCLUSION 

With the continuing importance of solar research, the use of libration point 
orbits between the Sun and the Earth is both an interesting and valuable area of 
effort. The need for orbit determination error analysis in conjunction with 
pre-mission station-keeping simulations was the original driving force behind this 
work. The results of three error analysis methods were compared with other 
similar libration point studies. The outputs of the error analysis were the six 
states’ standard deviations. These error levels could then, in turn, be used as 
error sources in Monte Carlo simulations of derived station-keeping routines. 
With nominal paths that could be constructed as nearly periodic halo-type, or 
distinctly quasi-periodic and smaller Lissajous trajectories, the error analysis 
and station-keeping results may differ by the type of orbit selected. Statistical 
tests for the equality of the average 2-year control costs using halo-type and 
Lissajous paths strongly suggest that there is no difference in mean 
station-keeping costs. It should, however, be noted that the results are 
presented for only one particular control algorithm and for two specific nominal 
trajectories. Experimental design methods are then used to determine the 
approximate functional relationship between the input uncertainties and the output 
2-year control cost. This type of functional relationship seems more useful than 
a series of tabular entries of control costs, each corresponding to a different 
set of input error levels. 
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